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L  INTRODUCTION 

The  field  of  three-dimensional  (3D)  underwater  acoustics  has  been 
developed  in  an  attempt  to  describe  the  properties  and  characteristics  of 
sound  propagation  in  the  heterogeneous  ocean  environment  (Lynch  and 
Chiu,  1989).  At  present,  there  are  three  major  approaches  used  in  studying 
the  3D  underwater  sound  field.  One  of  these,  ray  theory,  is  premised  upon  a 
high-frequency  approximation  method.  This  method  gives  a  geometric- 
optics  solution  which  involves  ray-tracing  in  a  spatially  varying  sound  speed 
field.  The  Hamiltonian  Acoustic  Ray  Tracing  Program  for  the  Ocean 
(HARPO)  is  a  versatile  numerical  code  developed  by  Jones  et  al.  (1986)  for  the 
computation  of  3D  rays.  The  primary  deficiency  of  this  theory  is  that  it  cannot 
adequately  address  the  behavior  of  low-frequency  sound  due  to  the  neglect  of 
sound  dispersion  and  diffraction. 

The  second  approach  uses  a  parabolic  approximation  to  the  acoustic  wave 
equation,  which  was  introduced  by  Tappert  (1977).  A  3D  numerical  parabolic 
equation  (PE)  model  was  developed  by  Lee  et  al.  (1988)  using  an  implicit  finite 
difference  scheme.  Another  3D  PE  model  was  arrived  at  by  Baer  (1981) 
utilizing  a  split-step  Fourier  algorithm.  In  a  recent  study,  analytic  solutions  to 
the  3D  parabolic  equation  were  obtained  by  Seigmann  et  al.,  (1990).  These 
solutions  are  valuable  for  testing  the  accuracy  of  3D  numerical  models. 

A  third  approach  is  the  normal  mode  method.  Pierce  presented  a  three- 
dimensional  version  of  this  method  in  1965.  The  normal  mode  method  is 
based  upon  a  separable  solution  to  the  wave  equation.  Pierce  assumed  an 
adiabatic  acoustic  environment  which   leads   to  the   neglect  of  coupling 


presented  by  Chiu  and  Ehret  (1990).  This  later  3D  normal  mode  model 
accounts  for  both  horizontal  refraction  and  model  coupling. 

Any  improvement  in  acoustic  modeling  must  be  quantified.  We  shall 
endeavor  to  test  the  accuracy  of  the  3-D  coupled  normal  mode  model  of  Chiu 
and  Ehret  (1990)  in  regard  to  horizontal  refraction. 

To  examine  the  accuracy  of  the  Chiu-Ehret  model,  we  compare  the  results 
of  their  model  with  those  from  analytic  solutions  to  the  parabolic  equation 
arrived  at  by  Seigmann,  et  al.  (1990).  The  results  of  the  Chiu-Ehret  model  and 
the  analytic  solutions  of  the  parabolic  equation  are  compared  for  two  cases  of 
horizontal  variation  in  sound  speed.  As  a  means  of  comparing  the  results,  we 
examine  the  slow  variations  of  the  complex  pressure  envelope  function  and 
transmission  loss. 

Our  simulated  acoustic  field  has  a  range  of  100  km  with  azimuth  from  30° 
to  180°.  Two  variations  in  the  horizontal  sound  speed  field  were  used.  In 
Case  I  the  sound  speed  varies  only  with  the  azimuth  angle  while  Case  II  varies 
both  with  the  azimuth  and  radial. 

The  normal  mode  theory  used  by  Chiu  and  Ehret  is  described  in  detail  in 
the  Appendix.  The  appendix  also  includes  a  description  of  the  PE 
approximation  and  the  analytic  solutions  developed  by  Seigmann  et  al. 
(1990). 


II. 


ANALYSIS  AND  METHOD  OF  MODEL  COMPARISON 


The  three-dimensional  coupled  normal  mode  model  was  applied  to  two 
sound  speed  fields  for  which  exact  analytic  solutions  to  the  PE  approximation 
are  available.  These  analytic  PE  solutions  for  vertically  invariant  sound  speed 
were  developed  by  Seigmann  et  al.  (1990).  Their  development  was  based  on 
expressing  the  "time-independent"  acoustic  pressure  component  as 
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C(f)  =  sound  speed  as  a  function  of  position   f. 

Substitution  and  separation  of  real  and  imaginary  parts  gives  two  coupled 
equations  governing  A   and  0: 
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If  one  chooses  @{r,6)  =  a(r)6,  equation  (2.2)  may  be  solved  for  the  amplitude, 


i.e., 


\(r,e)  =  F 


1  rra(s) 


ds 


(2.4) 


Choosing  a(r)  and  the  functional  dependence  F  results  in  a  reduced 
envelope  function  (and  ultimately  acoustic  pressure)  while  equation  (2.3)  can 
be  solved  for  the  index  of  refraction  n(f)  and  thus  the  sound  speed  field  C{f). 
Equation  (2.3)  is  for  a  wide  angle  PE  approximation.  This  wide  angle  form 
will  be  used  throughout.  To  judge  the  accuracy  of  the  Chiu-Ehret  model  in 
calculating  horizontal  refraction  we  compare  the  numerical  coupled  mode 
results  with  the  analytic  PE  results  for  two  sound  speed  fields  or  two 
functional  definitions  of  <x{r). 


A.    INDEX  OF  REFRACTION 

1.     Case  I 

For  Case  1  we  choose       a(r)  -  aQk0  r,  or 


0  =  a0k0  rO 


(2.5) 


where  aQ  is  an  arbitrary  constant.   From  equation  (2.4)  the  amplitude  is 

Ap(r,0)  =  F(0-aolnr). 

Let  the  function  F  be  multiplied  by  PQk0.  The  amplitude  has  the  final  form 

Ap(r,6)  =  P0k0(6-a0]nr)  (2.6) 

where  /3  0  is  an  arbitrary  constant  with  units  of  length.  The  square  of  the  index 
of  refraction  from  (2.5)  and  (2.6)  is 
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The  resultant  sound  speed  is  shown  in  Figure  1  for  a0  =  0.005  and 
C0  =  1500  m/s. 

2.     Case  II 

For  Case  II  we  choose 
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and  hence 
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The  square  of  the  index  of  refraction  is  now 
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The  resulting  sound  speed  is  shown  in  Figure  2  for  a0  =  4x10     and  C0  =  1500 
m/s. 


B.     ENVELOPE  FUNCTION 

We  now  have  two  analytic  sound  speed  fields  and  the  associated  acoustic 
pressure  fields  derived  from  the  parabolic  equation  approximation. 
Considering  only  j  =  1,  Seigmann  et  al.,  (1990)  define  pressure  as 
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For  Case  I,  using  the  amplitude  and  phase  expressed  in  (2.6)  and  (2.5), 
respectively,  we  obtain 
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For  Case  II,  using  (2.8)  and  (2.9),  we  get 


Figure  1.  Case  I  Sound  Speed  Field  for  a0  =  0.005  and  C0  =  1500  m/s 
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Figure  2.  Case  II  Sound  Speed  Field  for  aQ  =  4x10     and  C0  =  1500  m/s 
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In  terms  of  normal  modes,  pressure  can  be  expressed  as 


<PNM(r,e,z)  =  ^rY,Un(r,e)eil°K(r'e)d'Zn(z:r,e)  (2.14) 

th 

where  Zn  is  the  w  '  normal  mode  (eigenmode),  kn  is  the  corresponding 
horizontal  wavenumber  (eigenvalue)  and  Un  is  the  corresponding  slowly 
varying  modulation  envelope.   See  the  appendix  for  details. 

Since  we  are  considering  depth-invariant  sound  speed  fields,  an  analytic 
solution  is  available  for  the  eigenvalues  and  eigenmodes.    The  first  mode  is 

zi=  psin(7'2)  where  fi  =  V*2-*i  =^77 


and  thus,  keeping  only  the  contribution  of  the  first  mode,  (2.14)  becomes 

^NM^e-z)  =  jI-lu1(r/ey'o^(r-fl)''rsin(riz).  (2.15) 

The  normal  mode  model  entails  a  numerical  calculation  of  the  slowly 
varying  envelope  function  U  .  Therefore,  the  quantification  of  model 
accuracy  can  be  achieved  by  comparing  the  numerical  Un  to  those  derived 
from  the  analytic  PE  solutions.  The  analytic  envelope  function  11^,6)  for 
Case  I  is  derived  by  equating  (2.12)  to  (2.15).  The  resulting  expression  is 
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With  the  index  of  refraction  n  given  by  (2.7),  the  range  integral  of  /:-,  can  be 
expressed  as 


jk^dr  =kor< 


1  +  2 


i+ifn] 


-J4 


->2 


2 


-2 


4«0^  +  2ao  + 


^0  7 


It  follows  that  the  analytic  envelope  function  UA(r,0)  can  be  recast  as 


H 


Ul(r/e)=l—p0k0(e-a0\nr)c 


k0r<t>~ 


(2.17) 


where 


0  =  1 


4Uoy 


+  «o#-< 


3-J4 


1  +  I 
2 


Vk0J 


-2 


4a06  +  2aQ  + 


Kk0J 


For  Case  II,  the  wave  number  integral,  using  (2.10),  becomes 


jhdr  =  k0r 


2     2 


2 


l 

+— 

4 


+  a0k0re  +  -a$k$r2 


The  analytic  envelope  function  L71  is  then 


U1(r,e)  =  J-^(30k0(d-a0k0r)e 


k0r<p~ 


(2.18) 


where 


0=1- 


1-1 

4 


+  ccnknrO [  k-[dr 

VJ0 


For  Case  I  we  use  the  analytic  envelope  function  (2.17)  to  define 
the  initial  condition  at  r  =  r0  for  the  numerical  model  run.    The  condition  is 


Ui{r0,e)  =  J—pok0(0-a0\nr0)e^  4 


n 


(2.19) 


where 


0=i-_ 


irr,A 


v^oy 


l- 


+ao0- 


3-J4 


1  +  - 

2 


1^ 


-2 


4a06  +  2aQ  + 
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For  Case  II  we  use  the  analytic  envelope  function  from  (2.18)  at  r  =  rQ.  The 
Case  II  initial  condition  is 

Ul(r0,0)  =  J—  poko(0 -  aQk0r0)e  ^  (2.20) 


where 


<P  =  —7a0k0r0' 
6 
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III.  MODEL  RUNS 

Our  computational  domain  is  in  a  cylindrical  coordinate  system.  The 
horizontal  range  is  from  1  km  to  100  km.  The  azimuth  range  is  30°  to  150°. 
The  depth  is  4000  m.   The  domain  is  the  same  for  both  cases. 

We  use  a  frequency  of  50  Hz  and  source  depth  of  1000  m.  The  reference 
sound  speed  is  1500  m/s.  We  consider  mode  1  only.  The  vertical  boundary 
conditions  for  the  problem  are  pressure  release  surface  and  rigid  bottom.  Two 
range  and  azimuth  dependent  sound  speed  fields  as  presented  in  the  last 
chapter  are  used.  For  Case  I  sound  speed  ranges  from  1481  m/s  to  1496  m/s, 
varying  only  in  azimuth  with  a  gradient  of  0.128  to  0.130  m/s/degree  (see 
Table  1).  For  Case  II  sound  speed  varies  in  both  range  and  azimuth  (see  Table 
1  again).  Here  sound  speed  goes  from  1437  m/s  to  1500  m/s  with  a  radial 
gradient  of  0.13  to  0.657  m/s/km  and  an  azimuth  gradient  of  0.004  to  0.431 
m/s/degree.  Both  selections  of  sound  speed  fields  contain  realistic  gradients 
as  observed  in  the  ocean. 


TABLE  1. 

THE  AZIMUTH  AND  RADIAL  SOUND  SPEED  VARIATION 

Case 

Range 
r[km] 

Sound  Speed  [m/s] 

Azimuth  Variation 

Radial  Variation 

fP/km] 
dr    s 

6=30° 

6=90° 

9=150° 

6=30° 

6=90° 

e=i5o° 

0=30° 

0=90° 

0=150° 

I 

1 

1496 

14SS 

1481 

0.1303 

0.1289 

0.1275 

0 

0 

0 

50 

1496 

1488 

1481 

0.1303 

0.1289 

0.1275 

0 

0 

0 

100 

1496 

1488 

1481 

0.1303 

0.1289 

0.1275 

0 

0 

0 

II 

1 

1500 

1500 

1499 

0.0044 

0.0043 

0.0043 

0.1316 

0.3944 

0.6482 

50 

1493 

1480 

1468 

0.2217 

0.2178 

0.2142 

0.1310 

0.3848 

0.6221 

100 

1487 

1462 

1437 

0.4265 

0.4127 

0.4001 

0.1301 

0.3763 

0.5992 
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Eigenmodes  and  eigenvalues  are  calculated  using  finite  difference 
approximations  to  (A. 8)  in  the  appendix  with  a  40  m  vertical  grid  spacing  and 
a  matrix  eigenvalue  solver. 

To  obtain  the  numerical  envelope  function  a  first  order  differential 
equation  ((A. 24)  in  the  appendix)  is  integrated.  The  radial  integration  step 
size  is  0.5  km.  The  61  integration  paths  were  separated  by  two  degrees. 


13 


IV.    RESULTS 


A.    CASE  I 

In  comparing  the  numerical  normal  mode  (NM)  and  analytic  parabolic 
equation  (PE)  results  we  focus  on  the  slowly  varying  modulating  pressure 
envelope.  The  removal  of  the  rapid  oscillations,  which  have  wave  lengths  of 
order  27r/k0,  makes  the  displays  of  results  easier  to  interpret. 

Figures  3  and  4  show  the  amplitude  and  phase,  respectively,  of  the 
envelope  function  calculated  by  the  normal  mode  model  using  the  Case  I 
sound  speed  field.  Figures  5  and  6  are  the  amplitude  and  phase,  respectively, 
of  the  envelope  based  on  the  analytic  solution  to  the  parabolic  equation. 
Figure  7  shows  the  relative  difference,  defined  as 

UpE(  analytic  solution)  -  UNM  (numerical  solution) 
UpE  (analytic  solution) 

The  percent  difference  is  everywhere  less  than  5%.  The  difference  in  NM  and 
PE  phase  is  shown  in  Figure  8. 

Transmission  loss  at  a  depth  of  1000  m  calculated  by  the  NM  model, 
TLNM,  is  shown  in  Figure  9.  The  magnitude  of  the  analytic  PE  pressure,  from 
(2.12),  is 


|P|  =  J5VmM 


(3.1) 


Substituting  (2.6)  in  (3.1),  we  get 

\P\  =  J^-P(M0-oc0\nr)sm 
V  nknr 


n 
2H 
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Figure  3.   Amplitude  of  UNM  for  Case  I  with  a  =  0.005  and  0  =  0.1 


Figure  4.  Phase  of  UNM  for  Case  I  with  a  =  0.005  and  /3  =  0.1 
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Figure  5.   Amplitude  of  UPE  for  Case  I  with  a  =  0.005  and  ft  =  0.1 


Figure  6.  Phase  of  UpE  for  Case  I  with  a  =  0.005  and  0  =  0.1 
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Figure  7.  Amplitude  Difference  in  Percent  for  Case  I  with  a  =  0.005  and 

0=0.1 


Figure  8.  Phase  Difference  in  Degree  for  Case  I  with  a  =  0.005  and  ft  =  0.1 
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The  analytic  PE  transmission  loss  TLpE  is  thus 
TLPE  =  -201og10|P| 

=  101og10^-201og10|/3o/co(^  -  «0lnr)|  -  201og 


71 

sin z 

2H 


Figure  10  displays  TLpE  at  a  depth  of  1000m.  Both  transmission  loss  results 
show  that  the  azimuth  variation  of  sound  speed  produces  lower  loss  at  large 
azimuth  angle. 

The  difference  between  the  analytic  PE  transmission  loss  TLpE  and  the 
numerical  NM  transmission  loss  TLNM  is  shown  in  Figure  11.  Since  the 
modal  transmission  loss  is  a  function  of  amplitude  alone  this  difference  has 
the  same  shape  as  the  relative  error  in  the  amplitude  of  the  envelope  function. 
The  difference  in  transmission  loss  is  everywhere  less  than  1  dB. 

B.     CASE  II 

In  the  second  case,  we  examined  a  sound  speed  field  that  varies  in  range 
and  azimuth,  closely  representing  an  eddy  or  ring  structure  in  the  real  ocean. 

Figures  12  and  13  show  the  amplitude  and  phase,  respectively,  of  the 
envelope  function  UNN1  from  the  normal  mode  model.  Note  that  azimuth 
variation  of  the  amplitude  is  larger  than  range  variation.  Figures  14  and  15 
are  the  amplitude  and  phase  of  the  envelope  function  UpE  based  on  the 
analytic  solution  to  the  wide-angle  PE.  Figure  16  shows  the  relative 
difference  of  the  amplitude.  The  percent  difference  is  everywhere  less  than 
2%.  The  difference  between  the  NM  and  PE  phases  is  shown  in  Figure  17. 
The  difference  of  phase  is  everywhere  less  than  one  degree. 
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Figure  9.  TLNM  at  a  Depth  of  1000  m  for  Case  I  with  a  =  0.005  and  p  =  0.1 


Figure  10.  TLpE  at  a  Depth  of  1000  m  for  Case  I  with  a  =  0.005  and  p  =  0.1 
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Figure  11.  Difference  in  TL  in  dB  for  Case  I  with  a  =  0.005  and  j3  =  0.1 

The  NM  and  PE  transmission  losses,  TLNM  and  TLpE,  at  a  depth  of  1000  m 
are  shown  in  Figures  18  and  19,  respectively.  They  are  nearly  the  same  as  can 
be  seen  in  Figure  20  where  the  error  is  everywhere  less  than  0.2  dB. 
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Figure  12.  Amplitude  of  UNM  for  Case  II  with  a  =  4x10    and  (5  =  0.1 


-7 


Figure  13.  Phase  of  UNM  for  Case  II  with  a  =  4x10    and  p  =  0.1 
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Figure  14.  Amplitude  of  UPE  for  Case  II  with  a  =  4x10    and  fi  =  0.1 


Figure  15.  Phase  of  UPE  for  Case  II  with  a  =  4x10    and  (3  =  0.1 
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Figure  16.  Amplitude  Difference  in  Percent  for  Case  II  with  a  =  4x10    and  fi 

0.1 


-7 


Figure  17.  Phase  Difference  in  Degree  for  Case  II  with  a  =  4x10    and  j5  =  0.1 
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Figure  18.  TLNM  at  a  Depth  of  1000  m  for  Case  II  with  a  =  4x10"  and  (5  =  0.1 


Figure  19.  TLpE  at  a  depth  of  1000  m  for  Case  II  with  a  =  4x10'  and  p  =  0.1 
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Figure  20.  Difference  in  TL  in  dB  for  Case  II  with  a  =  4x10"  and  [3  =  0.1 
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V.        DISCUSSION 

In  this  chapter  we  address  two  issues.  The  first  is  a  quantification  of 
horizontal  refraction.    The  second  is  an  assessment  of  model  accuracy. 

A.    HORIZONTAL  REFRACTION 

Critical  to  this  accuracy  test  is  the  use  of  sound  speed  fields  that  cause 
significant  horizontal  sound  refraction.  To  confirm  that  the  selected  sound 
speed  fields  do  cause  significant  horizontal  refraction,  we  compare  the 
corresponding  Nx2D  NM  solutions  (the  calculation  is  divided  into  N  vertical 
slices  and  a  2D  NM  model  is  used  in  each  radial  direction)  with  3D  NM 
solutions. 

The  envelope  functions  of  the  N  x  2D  solutions  have  been  calculated 
from  (A. 25)  in  the  appendix.  The  amplitude  of  the  NM  envelope  function  is 
the  same  for  both  the  Nx2D  and  3D  solutions  in  Case  I  and  Case  II.  However, 
the  phases  are  different.  Figure  21  shows  the  phase  of  the  Nx2D  envelope 
function  for  Case  I.  The  phase  is  essentially  constant  at  -45°.  In  Case  II,  the 
phase  of  the  envelope  function  of  the  Nx2D  solution  is  also  constant  and  its 
value  is  also  -45°. 

Now  let  us  compare  the  phases  calculated  by  the  Nx2D  method  with  the 
phases  of  the  3D  solutions  for  Case  I  and  Case  II,  as  displayed  in  Figures  4  and 
13,  respectively.  Both  cases  show  a  13°-15°  phase  difference  over  a  range  of 
100  km,  implying  that  the  azimuthal  sound  speed  variation  used  is  large 
enough  to  induce  sound  propagation  out  of  the  vertical  plane. 
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Figure  21.   Phase  of  UNM  for  Nx2D  for  Case  I  with  a  =  0.005  and  £  =  0.1 

B.     ACCURACY 

1.     Phase 

In  this  section  we  discuss  the  differences  between  the  analytic  PE  and 
numerical  NM  solutions.  In  Case  I,  the  phase  difference  between  the  analytic 
PE  and  numerical  NM  solutions  is  less  than  2.8  degree  everywhere.  The  Case 
II  phase  differences  between  the  numerical  NM  and  analytic  PE  solutions,  as 
shown  in  Figure  17,  is  everywhere  less  than  1  degree.  A  quantification  of  the 
phase  differences  for  the  two  runs  is  given  in  Table  2. 


TABLE  2. 

PHASE  DIFFERENCE 

Case 

PE 

Normal  Mode  NM 

Difference  1  PE-NM  1 

min 

max 

min 

max 

max 

Case  I 

-45.15 

-60.35 

-45.15 

-58.92 

2.8 

Case  II 

-45.00 

-59.04 

-45.00 

-58.94 

0.43 
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2.     Amplitude 

Results  are  shown  in  Table  3.  Envelope  amplitude  of  NM  and  PE 
generally  agree  well  with  the  largest  relative  difference  near  5%.  The  largest 
difference  occurs  in  a  region  that  has  the  largest  azimuth  gradient  in  sound 
speed. 

TABLE  3.  AMPLITUDE  DIFFERENCE 


Case 

Amplitude  of  NM 

Amplitude  of  PE 

Error  Value 
Difference 

Relative 
Difference  [%] 

min 

max 

min 

max 

max 

max 

Case  I 

0.798 

4.22 

0.817 

4.27 

0.038 

4.8 

Case  II 

0.851 

4.28 

0.841 

4.28 

0.076 

1.8 

3.     Transmission  Loss  (TL) 

The  difference  in  transmission  loss  between  the  analytic  PE  and 
numerical  NM  solutions  is  shown  in  Table  4.  The  difference  for  both  cases  is 
everywhere  less  than  half  a  decibel. 

TABLE  4.  TRANSMISSION  LOSS  DIFFERENCE 


Case 

TL  of  NM 

TL  of  PE 

Difference 

min 

max 

min 

max 

max 

Case  I 

28.84 

63.31 

28.85 

63.73 

0.43 

Case  II 

28.73 

62.75 

28.73 

62.85 

0.15 
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VI.   CONCLUSIONS 

In  this  work  we  have  tested  a  three-dimensional  numerical  acoustic 
coupled  mode  model.  We  have  sought  to  describe  its  accuracy  in  regard  to 
the  modeling  of  horizontal  refraction.  In  order  to  test  its  accuracy,  we 
compared  the  coupled  mode  model  results  with  two  different  analytic 
solutions  to  the  parabolic  wave  equation. 

We  have  concentrated  on  the  accuracy  of  the  envelope  function  and 
transmission  loss  calculation.  For  the  sound  speed  fields  described  in  Table 
1,  we  found  that  the  phase  error  of  the  slowly  varying  envelope  function  is 
lower  than  2  degrees  and  that  the  envelope  amplitude  agrees  to  within  5%  or 
better.  Our  test  cases  have  indicated  that  the  normal  mode  model  agrees 
closely  with  the  analytic  solutions. 

We  have  only  used  depth-invariant  sound  speed  fields  for  this  test.  As  a 
result,  the  accuracy  in  modeling  mode-mode  interactions  is  not  tested  here. 
Future  tests  should  be  directed  at  examining  the  accuracy  of  calculating 
mode-coupling  effects  using  depth  varying  sound  speed  fields. 
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APPENDIX 

A.    THREE-DIMENSIONAL  COUPLED  NORMAL  MODE  METHODS 
1.     Overview 

a.      Three-Dimensional  (3D)  Helmholtz  Equation. 

The  acoustic  pressure  P(f,  t)  in  the  ocean  is  governed  by  the  wave 
equation: 


1     FP(T,t) 

C2{r)      dt 


V2P(r,t)-     ;"  'V     -0  (A-1) 


where  C(f)  represents  the  speed  of  sound  propagation.     The  cylindrical 
coordinates  are     f  =  (r,  6,  z)  where  r  is  range,  6  is   the   azimuthal   angle 

2 

measured   positive   counterclockwise,  -z  is  depth,  and  V     is  the  Laplacian 
operator  in  cylindrical  coordinates, 

^2    a2    i  a     i  a2     a2  ,A_. 

v  V  +  ^  +  ^?  +  ^  (A-2) 

Inclusion  of  an  acoustic  source  q{r,t),  modifies  Equation  (A.l)  to 


V2p(f,t)-±?-E^l  =  -4xq(r,t).  (A3) 


cz    a/z 


For  an  acoustic  point  source  located  at  r0   with   time-harmonic   circular 
frequency  w,  Equation  (A. 3)  becomes 

dt: 
If  we  let 


V2p(rj)--^d  ^  =  -44-f0)e-ia".  (A4) 
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p[r,t)  =  0(fyiwt  (a.5) 

where  Cf>(r)  is  time-independent,  then  a  Helmholtz  equation  results: 

V20(r)  +  k2{f)0{r)  =  -An8(f  -  f0)  ( A.6) 

w 

where  the  acoustic  wavenumber  is  k{f)  =  ——  .    For  r*  f n,  the  Helmholtz 

c(r)  ° 

Equation  is 

V2&(f)  +  k2(f)&{r)  =  0.  (A.6a) 

b.      Coupled  Normal-  Mode  Solution 

Adiabatic  normal  mode  theory  was  applied  to  a  range-dependent 
medium  by  Pierce  [1965]  and  Milder  [1969].  The  pressure  field  was  expressed 
as  a  linear  combination  of  the  local  modes  whose  coefficients  were  obtained 
from  a  system  of  decoupled  ordinary  differential  equations.  The  physics  of 
coupled  normal  mode  theory  includes  non-adiabaticity  and  results  in  a 
coupled  system  of  differential  equations.  The  theory  used  in  the 
development  of  the  Chiu-Ehret  3D  coupled  normal  mode  model  follows. 

The  first  step  in  the  development  is  the  expansion  of  the  time- 
independent  pressure  component  in  terms  of  the  local  normal  modes,  i.e., 


®(?)  =  ^yLPn(r,e)Zn(z:r,e)  (A.7) 


where 


Zn    is  the  n     local  normal  modes  at  point  (r,6). 

L  u 

Pn    is  the  n     mode  amplitude  function  at  point  (r,6). 
c.      The  Local  Normal  Modes  Zn 

The  local  normal  modes  obey  the  following  equation: 
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dr 


+ 


{k2-k2n) 


z„=o 


(A.8) 


th 


where  kn  is  the  horizontal  wavenumber  associated  with  the  n     mode  at  each 
horizontal  location. 

The  idealized  upper  and  lower  boundary  conditions  are 

Zn(z  =  O;r,0)  =  O.  (A.9) 


—  ZM(z  =  -H;r,0)  =  O. 

oz 


(A.10) 


where  H  is  the  ocean  depth.    The  local  normal  modes  can  be  normalized 
according  to  the  orthonormal  condition 

H 


J  ^n^m"2-  ~  °mn- 


(A.ll) 


0 


d.     The  Mode  Amplitude  Function  Pn 

To  obtain  the  governing  equation  for  P  ,  substitution  of  (A. 7)  in 
(A. 6a)  is  necessary.   The  use  of  (A.8)  and  the  farfield  approximation 


v1      11 


'  m^m  ~  *m'  m^m 


(A.12) 


following  the  substitution  results  in  the  following  equation: 


I 


m 


(d2P  1   d2P    ^ 

u  1  m   i     L   u  1  m 


■m 


dr 


+ 


dp 


+  2 


dP„dZ„       1  dPmdZ„\ 


m  "^m 


m  ^^m 


dr     dr       r2   dd    dd  ) 


+P, 


o  7Lm       1  o  Z 


w 


9r< 


+ 


m 


r2   dd2 


+  ^m^tn^m  ~  0- 


(A.13) 


The  next  step  is  to  multiply  (A.13)  by  Zn  and  then  integrate  over  z,  i.e., 
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I 


m 


jz„zm  L^L+^ii+4pmU+2jzi 

I    dr         r     dv~  I  V 


dPm  dZ„      1  3P,„  3Z 


m  ^^m 


+ 


m  ^^m 


dr     dr       r2   d0    dO  J 


dz 


J  *  m^n 


rd2Zm       1  d2Z 


+ 


m 


dz 
dr2       r2   dO2 


=  0. 


(A.14) 


A  subsequent  application  of  the  orthonormality  condition  (A. 11)  on  (A.14) 
gives 


,2      1    5 


2  \ 


(d2 


dr2   ■  "»      _2  ,*2 


r29^ 


=    I 


w 


2(JZ« 


azm  ,  >ap. 


jm 


3r 


dz 


n? 


+ 
dr 


4JZ« 


az 


m 


90 


dz 


3P. 


w 
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+ 


PmjZnl 


d  Zm  +l  d  Zm 


dr2       r2    dO2 


dz 


Defining  coupling  coefficients  as 


(A.15) 


y  =l\Z     dZm  (JZ 

J         dr 


dZ 


Pmn  ~        ^n 

r J  dO 


dz 


(d2Z 


\ 


dr 


1    7)27     ^ 


de2 


dz, 


the  coupled  system  (A.15)  can  be  recast  as 

^2  N 


(d2 


dr 


,2     J_  & 
'       n  +  r2d82 


Pn-1 


m  u 


1    d  \ 


Ymn  -n    +  Pmn     -,„ 
dr  r  dO 


+  B„,„  R 


mn 


m- 


(A16) 


e.     The  Envelope  Functions  Un 

Following  the  work  of  Chiu  and  Ehret  (1990),  the  mode 
amplitude  function  Pn  can  be  separated  into  a  slowly  varying  envelope 
function  Un  and  a  rapidly  varying  component  e     r'   ,  i.e., 


33 


Pn=Un{r,e)ei<p^'9) 

<Pn=\yn(r,0)dr. 
Jo 

A  substitution  of  (A. 17)  into  (A. 16)  yields 


(A17) 
(A18) 


rd2    ,2    i  a2 


dr2       "     r2d82 


Une^  =  -X 


m 


(      a    .    i d\ 


mn 


U  ei<Pm 


which  can  be  expanded  as 


fd2     i  a2^ 
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Dividing  (A.20)  by   i2kne  "  gives  the  governing  equation  for  the  envelope: 
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where 
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and  f  and  6  are  the  unit  directional  vectors  in  r  and  6,  respectively. 
/.      Transmission  Loss  TL 

The  mode  amplitude  function  P  n  is  related  to  the  envelope 
function  Un  and  the  phase  0M  by  (A. 17).  The  local  normal  modes  Zn  are 
solutions  to  the  eigenvalue-eigenfunction  problem  (A.8)-(A.ll).  Once  Pn  and 
Zn  are  computed,  pressure  can  be  calculated  from  the  product  of  Pn  and  Zn. 

.2 


At  each  point  (r,  6,  z),  <P  is 

,1,     -     *      r«    /*m2 


where 
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(A.22) 
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In  (A.22),  Re  and  Im  are  used  to  denote  the  real  and  imaginary 
parts,  respectively.   Transmission  loss  can  be  computed  as 


TL  =  -lOloglo<Z>O,0,z). 


(A.23) 
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2.      Numerical  Solution 

a.  Numerical  Model 

To  obtain  a  normal  mode  solution  we  must  solve  the  equations 
for  the  eigenvalues  and  eigenmodes  (A.8)-(A.ll)  and  the  equation  for  the 
modulation  envelope  (A. 21).   These  equations  can  be  solved  numerically. 

b.  The  Eigenvalues  k    and  the  Eigenmodes  Z|( 

Equation  (A. 8)  is  approximated  using  central  finite  differencing. 
This  approximation  casts  (A.8)-(A.10)  into  a  matrix  algebra  problem.  The 
eigenvalues  and  the  eigenmodes  of  which  can  be  determined  using  the 
iterative  QR  method  (Acton,  1970). 

c.  The  Iterative  Method 

To  obtain  the  envelope  function  U n,  we  need  to  solve  equation 
(A. 21)  which  can  be  rearranged  to  form  a  set  of  first  order  partial  differential 
equations  (PDE)  with  smaller  terms  put  on  the  right-hand  side.  We  can  solve 
for  U    in  an  iterative  fashion.  The  iterative  equation  is 
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where  Un  is  the  solution  at  the  /  iteration  for  mode  n.  This  set  of  first  order 
differential  equations  is  integrated  using  a  Runge-Kutta  method  of  order  5 
and  6  (Acton,  1970). 

3.     The  JVx  2D  Method 

In  the  N  x  2D  method,  the  sound  channel  is  divided  into  N  vertical 
slices  and  the  two-dimensional  normal  mode  solution  is  solved  for  each  slice. 
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It  combines  the  results  in  each  vertical  plane  to  construct  a  3D  field.    This 
method  provides  an  approximation  to  a  3D  solution  for  a  3D  sound  field. 

Thus,  in  each  slice  the  envelope  function  Un  is  independent  of 
azimuthal  variation  effects.  With  the  azimuthal  terms  ignored,  the  governing 
equation  becomes 
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B.     THREE-DIMENSIONAL  PARABOLIC  WAVE  EQUATION 
1.      Parabolic  Wave  Equation  (PE) 

The  following  work  was  developed  by  Tappert  (1977)  and  Lee  (1988). 
In  cylindrical  coordinates,  the  Helmholtz  Equation  (A. 6a)  becomes 


d20     1 d®      1  d2®     d2d>     ,2   2 
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(A26fl) 


where 
k0  =  w/C0 
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C0  is  a  reference  sound  speed 

w  is  the  harmonic  source  frequency 

n  is  the  index  of  refraction 
n  =  n(r,d,z)  =  C0  /C(r,9,  z). 

The  PE  approximation  assumes 

0(r,6,z)  =u(r,G,z)iir) 

where  v{r)  is  rapidly  varying  and  u(r,  0,  z)  is  a  modulation. 
By  substituting  (A. 26b)  into  (A. 26a),  one  obtains 
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If  v  is  determined  by 
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then  a  must  satisfy 
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The  solution  to  (A. 27)  is 

v(r)  =  til{kQr) 


ttVj 


exp 


for  -  — 
4 


tQ1 


(A.29) 


where  Ho  is  a  Hankel  function  of  zeroth  order. 

By  substituting  (A.29)  into  (A.28),  one  gets 
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which  can  be  factored  into 
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where 
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Q  is  called  the  square  root  operator.   If  one  considers  only  outgoing  waves,  the 
equation  to  be  solved  is 
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so  that 
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and  then  expand  Q  in  a  Taylor  series: 
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Inclusion  of  the  first  four  terms  in  the  expansion  leads  to  a  "wide  angle' 
version  of  the  PE  approximation.  Substituting  (A. 34)  into  (A. 32)  gives 
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2.     The  Envelope  Function  U    of  the  Parabolic  Equation 

Here  we  present  the  analytical  results  of  Seigmann  et  al.,  (1990).  Their 
development  uses  the  pressure  release  surface  and  rigid  bottom  boundary 
conditions,  i.e., 
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A  solution  to  equation  (A.35)  for  a  depth-independent  sound  speed  field  and 
subject  to  the  boundary  conditions  stated  above  can  take  the  following  form: 


u(r,e,z)  =  U;)(r,0)sin(y;z)exp 


ik0r 

2 

1 

1  -  - 

21 

2 

l*0j 

4 

[hi 

(A.37) 


where 


Yi 


.     1 


K 

77' 


H  is  the  ocean  depth,  and  U    is  the  envelope  function  of  the  parabolic 
solution. 

By  substituting  (A.37)  into  (A.35),  one  obtains  the  equation  governing 
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where  both  A  Jr,6)  and  &{r,Q)  are  real  quantities,  (A. 38)  becomes 

r  — 

2 


e^ 


ar 


a©      **(A> 


n2  - 1) 


1  + 


1 


yj 


K 


Oj 


Xf   2      ^2      1 
-  ?z     -1     -  - 

A  V  /  0 


2/cnr< 


2 

al4ii+    a^ae+      a^e    faef 
ae2  '  z  ae  ae  +2  v  oe2  "Ue. 


1*0  J 


40 


A  separation  of  the  real  and  imaginary  parts  of  (A.40)  yields 
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(A.42) 


Given  A  and  0,  a  pressure  field  is  determined,  and  the  index  of  refraction  can 
be  calculated  from  (A.42).  The  resultant  sound  speed  field  can  then  be  used 
in  the  normal  mode  model. 
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